Εργαστηριακή Άσκηση 2#
Σκοπός της δεύτερης σειράς ασκήσεων είναι η εξοικείωση με τις συναρτήσεις σχεδιασμού φίλτρων πεπερασμένης κρουστικής απόκρισης (FIR) και την υλοποίησή τους στο MATLAB. Προτού ξεκινήσετε την άσκηση θα πρέπει να μελετήσετε με προσοχή το Κεφάλαιο 1 και, ειδικότερα, την παράγραφο 1.3 του τεύχους του μαθήματος1. Το MATLAB (www.mathworks.com) είναι ένα διαδραστικό εμπορικό πρόγραμμα (Windows, Linux, Unix) με το οποίο μπορείτε να κάνετε εύκολα αριθμητικές πράξεις με πίνακες. Μπορεί να το έχετε εγκατεστημένο τοπικά, στον προσωπικό σας υπολογιστή ή να εργάζεστε σε κάποιο Εργαστήριο Προσωπικών Υπολογιστών (ΕΠΥ) της Σχολής σας που διαθέτει το συγκεκριμένο λογισμικό.
Μέρος 1: Εισαγωγή#
Στο MATLAB οι συναρτήσεις fft και ifft υποθέτουν ζεύγος μετασχηματισμού Fourier x(t) και X(f) υπολογισμένων σε μη αρνητικά διαστήματα t=[0:N-1]ts και f=[0:N-1]fo. Όπως έχετε ήδη δει στην Εργαστηριακή άσκηση 1, το άνω μισό μέρος του διαστήματος συχνοτήτων αντιστοιχεί στις αρνητικές συχνότητες του σήματος, όταν υπολογίζουμε το X(f) με τη βοήθεια της συνάρτησης fft (για σήματα πραγματικών τιμών, αυτά τα δυο μισά είναι κατοπτρικά ως προς το μέσο του διαστήματος). Ακριβώς το ίδιο ισχύει και για το άνω μισό μέρος του χρονικού διαστήματος, όταν το σήμα x(t) προκύπτει από τον αντίστροφο μετασχηματισμό Fourier μέσω της ifft. Το MATLAB διαθέτει τη συνάρτηση fftshift για να ολισθήσει κυκλικά τις τιμές του σήματος ή του μετασχηματισμού Fourier, ώστε να αντιστοιχούν σε κεντραρισμένα στο μηδέν αμφίπλευρα διαστήματα, δηλαδή, στις χρονικές στιγμές tb=[-ceil((N-1)/2): floor((N-1)/2)]ts ή στις συχνότητες fb=[–ceil((N-1)/2): floor((N-1)/2)]fo. Με τον τρόπο αυτό μπορούμε να παράγουμε τα xb(t) και Xb(f) που αντιστοιχούν στην αμφίπλευρη αναπαράσταση του σήματος και του μετασχηματισμού Fourier. Για να κατανοήσετε τα ανωτέρω θεωρείστε το διάνυσμα [1 2 3 4] ως το αποτέλεσμα του FFT μήκους 4. Τότε, το πρώτο στοιχείο (1) είναι ο όρος dc, το τρίτο στοιχείο (3) είναι το σημείο στο μισό της συχνότητας δειγματοληψίας fs/2, που μπορεί να εκληφθεί ότι αντιστοιχεί είτε στην –fs/2 είτε στην +fs/2. Τα στοιχεία 2 και 4 αντιστοιχούν στις συχνότητες +fs/4 και –fs/4. Εφαρμόζοντας την fftshift, το στοιχείο 3 εμφανίζεται πρώτο, που σημαίνει ότι στο MATLAB αντιστοιχεί στην αρνητική συχνότητα –fs/2, το επόμενο στοιχείο 4 αντιστοιχεί στη συχνότητα –fs/4 ακολουθούμενο από το dc και τη συχνότητα +fs/4. Για ένα μετασχηματισμό περιττού μήκους, δεν υφίσταται σημείο για το ±fs/2. Έτσι για το διάνυσμα [1 2 3], η εφαρμογή της fftshift θα δώσει τα στοιχεία που αντιστοιχούν στις συχνότητες –fs/3, 0, +fs/3. Εκτός του ότι παράγουν εξόδους με τις αρνητικές συχνότητες ή χρόνους στο άνω μισό του διανύσματος, αμφότερες οι συναρτήσεις fft και ifft αναμένουν ως είσοδο διάνυσμα με την ίδια μορφή, αφού προφανώς ισχύουν οι ταυτότητες h == ifft(fft(h)) και H == fft(ifft(H)) Η πρώτη υποδεικνύει ότι η είσοδος της ifft πρέπει να είναι αντεστραμμένη, όπως την παράγει η fft, και η δεύτερη ότι η είσοδος της fft πρέπει να είναι αντεστραμμένη, όπως την παράγει η ifft. Στο επόμενο σχήμα βλέπετε παραστατικά ένα ημιτονικό σήμα που έχει πολλαπλασιασθεί με παράθυρο Blackman τόσο στην αμφίπλευρη, όσο και την μονόπλευρη αναπαράστασή του.
Οι αντίστοιχοι μετασχηματισμοί Fourier, σε αμφίπλευρη και μονόπλευρη αναπαράσταση, φαίνονται στο επόμενο σχήμα:
Όταν τα x(t) και X(f) παράγονται από το MATLAB, δεν χρειάζεται κάποια ιδιαίτερη προσοχή, πλην της κυκλικής ολίσθησης σε περίπτωση που θέλουμε π.χ. να σχεδιάσουμε το αμφίπλευρο φάσμα ή σήμα. Όταν όμως ένα εκ των x(t) ή X(f) ορίζεται από τον χρήστη απαιτείται περισσότερη προσοχή, διότι, συνήθως χρησιμοποιούνται τα αμφίπλευρα σήματα ή φάσματα. Μπορείτε να μεταβείτε από τη μία αναπαράσταση στην άλλη ως εξής: x = ifftshift(xb), X = fft(x), Xb = fftshift(X), εάν ξεκινάτε από αμφίπλευρο σήμα και θέλετε να καταλήξετε σε αμφίπλευρο φάσμα, και X = ifftshift(Xb), x = ifft(X), xb = fftshift(x), εάν ξεκινάτε από αμφίπλευρο φάσμα και θέλετε να καταλήξετε σε αμφίπλευρο σήμα, όπου η συνάρτηση ifftshift του MATLAB εκτελεί την αντίστροφη λειτουργία της fftshift. Όταν το N είναι άρτιο, οι fftshift και ifftshift δίνουν το ίδιο αποτέλεσμα. Όταν όμως το N είναι περιττό αυτό δεν ισχύει και χρειάζεται προσοχή στη χρήση τους. Στην πράξη, η προσεκτική εφαρμογή των ανωτέρω έχει σημασία όταν υπολογίζεται η φάση του φάσματος. Το πλάτος του φάσματος δεν επηρεάζεται από την κυκλική ολίσθηση των στοιχείων που προκαλούν οι fftshift και ifftshift (δείτε ιδιότητες DFT).
Εξάσκηση#
Δοκιμάστε στο παράθυρο εντολών τα ακόλουθα προκειμένου να εμπεδώσετε τη χρήση των συναρτήσεων fftshift και ifftshift.
X=[-2:2] fftshift(X) ifftshift(X) Y = fftshift(fftshift(X)); Z = ifftshift(fftshift(X)); isequal(X,Y) isequal(X,Z)
Ερώτηση 1: Ποιο εκ των διανυσμάτων Υ και Ζ ισούται με το X; Γράψτε την απάντησή σας σε ένα αρχείο κειμένου lab2_nnnnn.txt, όπου nnnnn τα πέντε τελευταία νούμερα του αριθμού μητρώου σας, χρησιμοποιώντας το Notepad από το μενού των Windows (Start → Programs → Accessories → Notepad) και αποθηκεύστε το στον φάκελο My Documents. Θα υποβάλετε το αρχείο αυτό ηλεκτρονικά στο τέλος, αφού απαντήσετε και τις επόμενες ερωτήσεις, οπότε μπορείτε να τα αφήσετε ανοικτό. Ερώτηση 2: Επαναλάβατε με X=[-1:2]. Τι παρατηρείτε; Γράψτε την απάντησή σας στο αρχείο κειμένου lab2_nnnnn.txt. Δοκιμάστε στο παράθυρο εντολών τα ακόλουθα δύο παραδείγματα για να εμπεδώσετε τη χρήση των συναρτήσεων fftshift και ifftshift σε συνδυασμό με τις fft και ifft.
close all; clear all;clc; xb=[1 2 3 4 5 4 3 2 1] % πραγματικό σήμα με άρτια συμμετρία figure; subplot (2,1,1); plot([-4:4],xb); ylabel(‘xb’); x=ifftshift(xb) % το σήμα με τις αρνητικές συνιστώσες στο άνω μέρος X=fft(x) % FFT Xb=fftshift(X) % το φάσμα με τη dc συνιστώσα στο κέντρο, πραγματικές
subplot (2,1,2); plot([-4:4],Xb);ylabel(‘Xb’); close all; clear all;clc; Xb=[0 0 1 1 1 1 1 0 0] % φάσμα βαθυπερατού σήματος με άρτια συμμετρία figure; subplot (2,1,1); plot([-4:4],Xb); ylabel(‘Xb’); X=ifftshift(Xb) % το φάσμα με τις αρνητικές συνιστώσες στο άνω μέρος x=ifft(X) % IFFT xb=fftshift(x) % πραγματικό σήμα με άρτια συμμετρία όπως αναμένεται subplot (2,1,2); plot([-4:4],xb); ylabel(‘xb’);
Ερώτηση 3: Τροποποιείστε το προηγούμενο παράδειγμα ώστε να ξεκινήσετε απευθείας με τον ορισμό του φάσματος του βαθυπερατού σήματος X όπως το αναμένει η ifft. Γράψτε την απάντησή σας στο αρχείο κειμένου lab2_nnnnn.txt.
Show code cell source
from scipy import signal
import scipy.io.wavfile
# Ανατρέξτε στην τεκμηρίωση της βιβλιοθήκης scipy.signal
# https://docs.scipy.org/doc/scipy/reference/signal.html
from scipy.fft import fft, fftfreq
from dash import Dash, dcc, html, Input, Output
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from numpy import random
import pandas
import plotly
import plotly.express as px
import plotly.graph_objs as go
import dash_bootstrap_components as dbc
import warnings
warnings.filterwarnings('ignore')
print("Libraries added successfully!")
Libraries added successfully!
Μέρος2: Σχεδιασμός και υλοποίηση φίλτρων#
Θα ασχοληθούμε με το Παράδειγμα 1.2 της παραγράφου 1.5 του τεύχους Μαθήματος. Το παράδειγμα αυτό παρουσιάζει δύο εναλλακτικές μεθόδους σχεδιασμού FIR φίλτρων: α) τη μέθοδο των παραθύρων και β) τη μέθοδο των ισοϋψών κυματώσεων τις οποίες εφαρμόζει στην περίπτωση βαθυπερατών φίλτρων.
Στο παράδειγμα, τα φίλτρα δοκιμάζονται σε ένα πραγματικό σήμα, s, το οποίο είναι αποθηκευμένο στο αρχείο sima.mat (binary αρχείο MATLAB). Πρόκειται για ένα σήμα sonar με φάσμα που εκτείνεται μέχρι περίπου τα 4 KHz και συχνότητα δειγματοληψίας Fs=8192 (είναι και αυτή αποθηκευμένη στο αρχείο sima.mat, μαζί με το σήμα).
Εδώ θα πειραματιστούμε με δύο σήματα: (i) το sonar του παραδείγματος, το οποίο εδώ διαβάζεται από ένα .txt αρχείο (έχει προέλθει με εξαγωγή του s από το MATLAB) και (ii) ένα σήμα μουσικής, το violin.wav (σήμα από μουσική βιολιού), το οποίο περιέχει υψηλότερες συχνότητες και έχει προέλθει με δειγματοληψία στα Fs_viol=44100 Hz.
Σήμα sonar#
Show code cell source
# Ανάγνωση δειγμάτων σήματος από txt file
with open('files/sima.txt') as f:
s = [float(x) for x in f]
s=np.array(s)
print('μέγεθος σήματος=', s.shape)
Fs=8192
μέγεθος σήματος= (6565,)
Στο πεδίο του χρόνου#
…
Show code cell source
t=np.arange(0,len(s))/Fs
# Time domain plot of x
fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=s, mode='lines', line=dict(color='#00CC96')))
fig.update_layout(title='Time domain plot of x', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
xaxis_title='t (sec)',
yaxis_title='Amplitude',
template='plotly_white',
xaxis=dict(range=[0, 0.8]),
yaxis=dict(range=[-0.05, 0.05]))
fig.show()
Ακούμε το σήμα#
Show code cell source
# Πρέπει να έχουμε εγκατεστημένη τη βιβλιοθήκη sounddevice
import sounddevice as sd
import ipywidgets as widgets
from IPython.display import display
def play_sound(b):
sd.play(20*s,Fs)
# Create a button widget
play_button = widgets.Button(description="Play Sound")
# Link the button to the function that plays the sound
play_button.on_click(play_sound)
# Display the button
display(play_button)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[4], line 2
1 # Πρέπει να έχουμε εγκατεστημένη τη βιβλιοθήκη sounddevice
----> 2 import sounddevice as sd
3 import ipywidgets as widgets
4 from IPython.display import display
ModuleNotFoundError: No module named 'sounddevice'